Natural human mobility patterns and spatial spread of infectious diseases 
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We investigate a model for spatial epidemics explicitly taking into account bi-directional move- 
ments between base and destination locations on individual mobility networks. We provide a system- 
atic analysis of generic dynamical features of the model on regular and complex metapopulation net- 
work topologies and show that significant dynamical differences exist to ordinary reaction-diffusion 
and effective force of infection models. On a lattice we calculate an expression for the velocity of the 
propagating epidemic front and find that in contrast to the diffusive systems, our model predicts a 
saturation of the velocity with increasing traveling rate. Furthermore, we show that a fully stochas- 
tic system exhibits a novel threshold for attack ratio of an outbreak absent in diffusion and force 
of infection models. These insights not only capture natural features of human mobility relevant 
for the geographical epidemic spread, they may serve as a starting point for modeling important 
dynamical processes in human and animal epidemiology, population ecology, biology and evolution. 
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The geographic spread of emergent infectious dis- 
eases, epitomized by the 2009 HlNl outbreak and subse- 
quent pandemic [Ij, the worldwide spread of SARS in 
2003 ^ and recurrent outbreaks of influenza epi- 
demics [IHSl, is determined by a combination of disease 
relevant human interactions and mobility across multiple 
spatial scales [7j. While infectious contacts yield local 
outbreaks and proliferation of a disease in single popula- 
tions, multiscale human mobility is responsible for spatial 
propagation [8j. Therefore, a prominent lineage of math- 
ematical models has evolved that is based on reaction- 
diffusion dynamics ^ [10] in which the combination of 
local exponential growth and diffusive dispersal captures 
qualitative aspects of observed dynamics. Typically these 
systems exhibit constant velocity epidemic wave fronts. A 
related class of phenomenological models is based on the 
concept of an effective force of infection across distance 
and thus does not require explicit modelling of disper- 
sal [TT|[T2]. 

Although more sophisticated models [U HI [131 (Hi have 
been developed to describe the dynamics of recent emer- 
gent epidemics such as the HlNl pandemic or SARS, tak- 
ing into account long distance travel and multiscale mo- 
bility networks [15j, the majority of these models are still 
based on the interplay of local reaction kinetics and diffu- 
sion processes on networks of metapopulations. The key 
assumptions of diffusive transport (Fig{l](a)) are that a) 
individuals behave identically, b) movements are stochas- 
tic, c) spatial increments are local and as consequence 
individuals eventually visit every location in the system. 
Although it is intuitively clear that these assumptions 
are idealizations in conflict with everyday experience, the 
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difficulty is how to refine them when data on mobility is 
lacking, is insufficient or incomplete. Fortunately, a se- 
ries of recent studies [TF-TO] substantially advanced our 
knowledge on multiscale human mobility. An important 
discovery that emerged from these studies, are individual 
mobility networks, i.e. individuals typically only visit a 
limited number of places frequently, predominantly per- 
forming commutes between home and work locations and 
possibly a few other locations. Consequently individu- 
als or groups of individuals exhibit spatially constrained 
movement patterns despite their potentially high mobil- 
ity rate. It has remained elusive how this novel and im- 
portant empirical insight on individual mobility networks 
can be reconciled with epidemiological models, to which 
extent it may impact spatial disease dynamics, and how it 
may alter spreading scenarios and predictions promoted 
by ordinary reaction-diffusion models, in which mobile 
hosts can reach every location in the system. 

In this Letter we demonstrate how individual mobil- 
ity networks can be incorporated into a class of models 
for spatial disease dynamics, and address the questions 
whether and how significantly the existence of spatially 
constrained individual mobility networks impact on key 
features of disease dynamics. To this end we investigate 
a model that explicitly accounts for bidirectional mobil- 
ity of individuals between their unique base location (e.g. 
their home) and a small set of other locations (Figjljb)). 
The entire population is therefore represented by a set 
of overlapping individual mobility networks associated 
with each base location, see e.g. [20-22J. We focus on 
the analysis of epidemics on regular lattices and com- 
plex metapopulation networks and systematically com- 
pare the dynamics to reaction-diffusion systems as well 
as force of infection models. In conflict with reaction- 
diffusion models, in which front velocities increase with 
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travel rates unboundedly, we show that our model pre- 
dicts a saturation of wave front velocities, a direct con- 
sequence of the rank of locations in individual mobility 
patterns. This suggests that estimates for propagation 
speeds may have been considerably overestimated in the 
past. Furthermore, we analyze a fully stochastic model 
to show that both, in regular lattices as well as complex 
metapopulation networks the global outbreak of a disease 
is determined by a novel type of threshold for the attack 
ratio of the disease which depends on the charactetistic 
time spent at distant locations. Finally, we show that in 
the limit of low and high travel rates our model agrees 
with reaction-diffusion models and the direct force of in- 
fection class of models, respectively. 

Consider a system of M populations labeled m and as- 
sume that in each an epidemic outbreak can be described 
by a compartmental SIR-model 

Im H~ ^ 2/7^, Iyyi y Rm-) (1) 

where 5^, and R^n label and quantify susceptible, in- 
fected and recovered individuals of population m. In- 
fections and recovery events occur at rates a and /3, re- 
spectively, with a > j3. The number of individuals in 
population m is given by N^^ = Sm + /m + Rmi which 
though is conserved only in statistical sense at equilib- 
rium. This system of reactions yields the mean-field de- 
scription dtl = alS/N - pi and dtS = -alS/N. A 
natural and plausible extension to a system of M cou- 
pled populations is diffusive dispersal among those pop- 
ulations defined by a hopping rates Wnm > from pop- 
ulation m to n, which yields a metapopulation reaction- 
diffusion system, i.e. (for infectives and susceptibles) 

m 

dtSn = -aSnln/N^ + ^ (^nm^'m - ^mn^n) (2) 
m 

where is the number of individuals in population n 
in diffusive equilibrium. Travel rates uornn are usually 
estimate by gravity-like laws [23j. is determined by 
detailed balance, i.e. N^/N^ = Wnm/'^mn and conserva- 
tion of the number of individuals in the metapopulation 
M = . These type of models have been em- 

ployed in numerous recent studies [2j [Ml [241 IIS]- The 
tight connection to spatially continuous reaction diffu- 
sion systems is revealed for the special case of a linear grid 
of populations placed at regular intervals / at positions 
Xm = ml and mobility between neighboring populations. 
I.e. WjijYi = ix)5n-i^m + ^^n+i,m which yields 

dtj^ajs-pj^Ddlj (3) 
dtS = —ajs + Dd^s 

in which j{x,t) = In/N^, s{x,t) = Sn/N^ and D = Pw. 
Eq. (|3| is related to a classic form of a Fisher equa- 
tion which has been deployed in mathematical epidemi- 
ology ^[lOj. For sufficiently localized initial conditions 




Figure 1. (color online) Metapopulational mobility mod- 
els. Patches and arrows represent individual populations and 
travel, (a) Diffusive dispersal: indistinguishable individuals 
travel randomly between different locations governed by the 
set of transition rates Unm- (b) Dispersal capturing individ- 
ual mobility patterns. Individuals with label k possess travel 
rates cj^/. and cj^^ at which they travel from their base loca- 
tion k to connected locations m and back. 

this systems exhibits travelling waves with speed 

c = 2^aD{l- p/a) - V^. (4) 

Note that this velocity increases monotonically with the 
mobility rate uo. 

In order to account for individual mobility patterns we 
propose the following generalization: individuals possess 
two indices n and /c, characterizing their current and their 
base location, respectively. The dispersal dynamics is 
governed by a Markov process 

21^ /K\ 

a; 

nm 

where X stands for each infectious state I and R. 
Eqs. ([5| imply that individuals with base location k pos- 
sess a unique mobility rate that determines how 
they travel between locations n and m. This yields a 
generalization of ([2|: 

9tln = l^^n ^^^n ~ P^n + X] i^nm^m ~ ^mn^n) 
m m 

^ m m 

where and are the number of infecteds and sus- 
ceptibles of type k that are currently located in popula- 
tion n. Nn is the total number of individuals at n, i.e. 
Nn = ^kiln^S^^R^). If uj^^ are /^-independent, we re- 
cover the reaction-diffusion case described above. In the 
following we consider the case of overlapping star-shaped 
networks corresponding to commuting between base and 
destination locations only. This implies that lo^^ = if 
either k ^ n or k ^ m. We further assume the constant 
return rate w^j^ = w~ for all k and m. Realistically we 
have uj]^^/uj~ <C 1, implying individuals belonging to n 
remain at their base most of the time. 

If mobility rates are large compared to the rates associ- 
ated with the infection and recovery, i.e. uj^^^uj~ ^ a, /3, 
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detailed balance is fulfilled for infecteds and susceptibles 
separately and the last term in Eq. (|6| vanishes in which 
case the above model is equivalent to the remote force of 
infection model, see also [11, 21j. 

In general, it is difficult if not impossible to extract 
the dynamic differences between the systems defined by 
Eqs. Q and (|6| for an arbitrary metapopulation. Yet, 
it is crucial to understand the key dynamic differences 
that emerge when individual mobility patterns replace 
ordinary diffusion. We therefore investigate the impact 
of individual mobility patterns in the instructive one- 
dimensional lattice case. We assume that only infecteds 
can travel and that recovery is absent (SI model), i.e. 
/3 = (relaxing these two restrictions does not change 
the main results but clarifies the analysis). We denote 
the number of infecteds with base location n that are 
located on site n (i.e. their base) and sites n ± 1 by 
and , respectively. At rates cj+ and uj~ indi- 
viduals leave and return to their base. In order to ob- 
tain a spatially continuous description we approximate 
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relative concentrations 
leading to 



^A/^ and introduce 
IJN, vn = {I++I-)/N 



dtu = a{l — u — v){u -\- V -\- VAv) -\- u v — 2uj~^u 
dtv = 2uj^u — uj~v (7) 

where V = /^/2. The function u{x^t) is the fraction of 
individuals based and located at x whereas v{x,t) is the 
fraction of individuals based but not located at x. Note 
that we discarded here a third equation for w = (/+ — 
I~)/N independent from ^ with a solution vanishing 
at long times. This does not change the main result 
The non-trivial steady state is given by u"^ = uo~ /(2uo^ 
cj~), v'^ = 2co'+/(2co'++co'~) with ^v'^ = 1 as expected. 

One of the key questions is if the system ^ ex- 
hibits stable propagating waves as solutions and how 
their velocity depends on system parameters. The ansatz 
uix^t) = u{x — ct) and likewise for leads to a stable 
solution with a velocity given by 
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(8) 



Fixing uj~ this implies that c ^ as uj^ 0. On the 
other hand if a;~ is small, the system is entirely deter- 
mined by forward rate uj^ . For a systematic comparison 
to the ordinary reaction-diffusion system, we have to es- 
tablish a relation between uo^ and uj in Eq. Q. For the 
sake of simplicity we consider the symmetric case, i.e. 
uj^ = uj~ = uj^ yielding 
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Numerical simulations of a fully stochastic system nicely 
agree with the analytical predictions (Fig.[2|. The essen- 
tial feature of c{uj) is its saturation as ^ oo whereas 
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Figure 2. (color online) Front velocity c{uj) of the model de- 
fined by Eq. (|6]) as a function of mobility rate uj in compar- 
ison to ordinary reaction-diffusion dynamics (Eq. For 
stochastic simulations we used a lattice with 10^ sites and 
N — lO'^ individuals/site. Mean velocities are indicated 
by red (bi-directional mobility) and blue (reaction-diffusion) 
symbols. Analytical results, i.e. Eqs. ^ and Q are shown 
by dashed and dash-dotted lines respectively. 



c ~ ^/uJ for reaction-diffusion systems. This effect is a 
direct consequence of the spatial constraints imposed by 
individual mobility patterns, i.e. increasing the mobil- 
ity rate uo yields a higher frequency of travel events but 
restricted to the set of locations connected to the base 
location and thus is universal and holds also for more 
complex metapopulation topologies [26j. For high rates 
the deviation between ordinary reaction diffusion and bi- 
directional mobility increases without bounds. This is 
an important insight as the expression for wave front 
speeds for reaction-diffusion have been used in the past 
to estimate wave front speeds from mobility rates and 
vice versa. Our results indicate that if constrained mo- 
bility patterns and bi-directional movement patterns are 
at work these estimates must be treated with particular 
care. Note that in the limit cj oo Eq. ^ reduces to 
the heuristic equation in Ref. [27j . 

A key question in epidemiological contexts concerns 
conditions under which an epidemic outbreak propa- 
gates or wanes. Outbreak dynamics is usually triggered 
by crossing thresholds inherent in the system's dynam- 
ics [28j. For example the basic reproduction number IZq 
of a disease, i.e. the expected number of secondary cases 
caused by a single infected individual in a susceptible 
population represents such a threshold [8j . In single pop- 
ulation SIR dynamics introduced above IZq = a/j3. If 
IZq > 1 an outbreak occurs, otherwise epidemic dies out. 
In the extended metapopulation system with diffusive 
host movements a second threshold, known as the global 
invasion threshold [24J , is induced by the flux of individ- 
uals between populations which must be sufficiently high 
for an epidemic to spread throughout the metapopula- 
tion. 

An interesting and important feature of the model we 
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Figure 3. Results of stochastic simulations. Attack ratio 
p{(jj~ /N) as a function of the return travel rate uj~ for SIR 
epidemics on a lattice with 10^ sites with N agents/site. Epi- 
demic parameters: a = 1, /3 = 0.1. Insets: (a) — p{uj~ /N) 
for a scale-free network (7 = 1.5) of 10^ nodes populated uni- 
formly with {N) — 250/site, with /cmin — 5 and /cmax — 50; 
(b) — p{uj~ /N) for an Erdos-Renyi network with 500 nodes 
and ^ = 10. Results were averaged over 50 realisations. In- 
terlocation flux rate was (cj) = 1. 



propose is the existence of another invasion threshold 
which is determined by the return rate or equivalently 
by the time an individual spends outside the base loca- 
tion. This finding is illustrated in Fig. |3] which depicts 
the attack ratio p as a function oi uj~ for a bidirectional 
SIR epidemic on a lattice. The attack ratio p is defined 
as the fraction of the overall population affected during 
an epidemic. We fixed the entire interlocation flux at a 
value above the global invasion threshold ensuring global 
outbreak in the reaction-diffusion system. 

This result is universal for bidirectional system and is 
independent of the particular topology. Insets (a) and (b) 
of Fig. [3] display simulation results for uncorrelated scale- 
free and Erdos-Renyi networks. For low return rates the 
attack ratio is close to unity. However, with growing val- 
ues of the return rate, the attack ratio decreases steeply 
and vanishes, reflecting the absence of the global out- 
break. The regime of high return rates corresponds to 
small dwelling times on distant locations. Infected in- 
dividuals do not have enough time to transfer the dis- 
ease to susceptibles in unaffected locations before they 
return. Note that for a wide range of local population 
sizes (50 — 2, 000) the attach ratio p collapses onto a uni- 
versal curve p ~ p{(jj~ /N) which suggests that the basic 
mechanism of the invasion threshold and its functional 
dependence on return rate uo~ can be understood theo- 
retically. 

To estimate the observed invasion threshold on a one- 
dimensional lattice for IZq ^ 1 we calculated the number 
of infecteds A that originate from an affected location and 
can seed an outbreak in an unaffected one [28]. These 



seeders are approximately given by the total number of 
individuals entering the new location per unit time, i.e. 
Nuo multiplied by the typical time they remain active in 
the destination location. This time is given by the inverse 
rate r of exiting the seeders class. Since two processes 
(returning to the base location and recovery) can trigger 
this exit, r = /3 + uj~ . Therefore A ~ Nuo/{j3 + cj~). 
In the tree approximation the threshold condition for an 
epidemic on a network with an average degree k is given 
by A(^ — l)(7^o — 1) > 1 [24J. Thus, for a one-dimensional 
lattice the threshold condition reads 



(10) 



where we kept the total interlocation flux uj = 
2uo^ uo~ / {2uo^ uj~) constant (this relation and positiv- 
ity of travel rates impose the restrictions: uo~ > uo and 
uj^ < uj/2) and used a;+ <C a;~ . Moreover, as extensive 
simulations show, the global invasion threshold in terms 
of the flux rate uj exists in bidirectional systems only for 
low return rates uj~ . With increasing return rate the 
disease fails to propagate globally, even for constant in- 
terlocation flux UJ. This effect, observed in lattice topolo- 
gies and paradigmatic more realistic network topologies 
is a fundamental consequence of birectional mobility pat- 
terns. This effect is absent in reaction-diffusion systems 
fully determined by the total interlocation flux. This 
property implies that invasion of a front is limited only 
by return rates and not as is typically assumed by the 
overall particle flux uo. In the context of disease dynam- 
ics this implies that more efficient containment strategies 
could be potentially devised that do not target the over- 
all mobility but rather modifies mobility patterns in an 
asymmetric way. 

Note that the transition takes place only in a system 
with a finite number of agents per site, in the system with 
an infinite number of individuals this effect would disap- 
pear as also confirmed in Fig. [sl Indeed from Eq. (10) 
for /3 <C cj" the scaling p ~ pt^^) follows — with in- 
creasing number of individuals per site the threshold 
value of the return rate uj~ increases. 

Recently an unprecedented amount of detailed infor- 
mation on human activity became available requiring re- 
vision of established models and new more sophisticated 
ones. We investigate an epidemiological model explicitly 
incorporating such important properties of human mo- 
bility as individual mobility networks and frequent bidi- 
rectional movements between home and destination loca- 
tions. It manifests surprising dynamical features such as 
the existence of bounded front velocity and novel propa- 
gation threshold which are universal for any metapopula- 
tion topology. As more detailed data continues to become 
available our approach is a promising foundation for the 
further research. 

D.B. and V.B. acknowledge support from the Volkswa- 
gen Foundation, D.B. acknowledges EU-FP7 grant Epi- 
work. V.B and D.B would like to thank W. Noyes, and 
H. Schlammer for invaluable comments. 



5 



[1] Christophe Fraser, Christl A Donnelly, Simon 
Cauchemez, et a/., "Influenza: Making privileged 
data public response," Science, 325, 1072-1073 (2009). 

[2] L Hufnagel, D Brockmann, and T Geisel, "Forecast and 
control of epidemics in a globalized world," Proc. Natl. 
Acad. Sci. USA, 101, 15124 (2004). 

[3] Vittoria Colizza, Alain Barrat, Marc Barthelemy, et aL, 
"Predictability and epidemic pathways in global out- 
breaks of infectious diseases: the sars case study," BMC 
Med., 5, 34 (2007). 

[4] D. Balcan, H Hu, B Goncalves, et a/., "Seasonal trans- 
mission potential and activity peaks of the new influenza 
A(HINI): a Monte Carlo likelihood analysis based on 
human mobility," BMC Med., 7, 45 (2009). 

[5] D Balcan, B Goncalves, H Hu, et a/., "Modeling the spa- 
tial spread of infectious diseases: The GLobal Epidemic 
and Mobility computational model," J. Comp. Sci., 1, 
132 (2010). 

[6] D Balcan, V Colizza, B Goncalves, et aL, "Multiscale 
mobility networks and the spatial spreading of infectious 
diseases," Proc. Nat. Acad. Sci. USA, 106, 21484 (2009). 

[7] S Riley, "Large-scale spatial-transmission models of in- 
fectious disease," Science, 316, 1298 (2007). 

[8] R M Anderson and R M May, Infectious Diseases of Hu- 
mans (Oxford University Press, 1991). 

[9] A Kolmogorov, I Petrovsky, and N Piscounov, "Etude 
de I'equation de diffusion avec croissance de la quantite 
de matiere et son application a un probleme biologique," 
Bull, de I'univ. d'etat a Moscou, Ser. internat., sect. A, 
1, 1 (1937). 

[10] R A Fisher, "The wave of advance of advantageous gene," 
Ann. Eugenics, 7, 355-369 (1937). 

[11] S Rushton and A J Mautner, "The deterministic model 
of a simple epidemic for more than one community," 
Biometrika, 42, 126-132 (1955). 

[12] T J Hagernaas, C A Donnelly, and N M Ferguson, "Spa- 
tial heterogeneity and the percictence of the infectious 
diseases," J. Theor. Biol., 229, 349 (2004). 

[13] L A Rvachev and I M Longini, "A mathematical model 
for the global spread of influenza," Math. Biosc, 75, 3 
(1985). 

[14] V. Collizza, R. Pastor- S at orr as, and A. Vespignani, 
"Reaction-diffusion processes and metapopulation mod- 
els in heterogeneous networks," Nat. Phys. (2007). 

[15] Alessandro Vespignani, "Predicting the behavior of 



techno-social systems," Science, 325, 425 (2009). 

[16] D Brockmann, L Hufnagel, and T Geisel, "The scaling 
laws of human travel," Nature (London), 439, 462 (2006). 

[17] Marta C Gonzalez, Cesar A Hidalgo, and Albert- 
Laszlo Barabasi, "Understanding individual human mo- 
bility patterns," Nature (London), 453, 779 (2008). 

[18] Chaoming Song, Zehui Qu, Nicholas Blumm, and 
Albert-Laszlo Barabasi, "Limits of predictability in hu- 
man mobility," Science, 327, 1018 (2010). 

[19] Chaoming Song, Tal Koren, Pu Wang, et a/., "Modelling 
the scaling properties of human mobility," Nat. Phys., 6, 
818 (2010). 

[20] Lisa Sattenspiel and Klaus Dietz, "A structured epidemic 
model incorporating geographic mobility among regions," 
Math. Biosc, 128, 71 (1995). 

[21] Matt J. Keeling and Pejman Rohani, "Estimating spa- 
tial coupling in epidemiological systems: a mechanistic 
approach," Ecol. Lett., 5, 20-29 (2005). 

[22] J Arino and P van den Driessche, "The basic repro- 
duction number in a multi-city compartmental epidemic 
model," Lect. Notes in Control and Information Science, 
294, 135 (2003). 

[23] V V Belik and D Brockmann, "Accelerating random 
walks by disorder," New J. Phys., 9, 54 (2007). 

[24] Vittoria Colizza and Alessandro Vespignani, "Invasion 
threshold in heterogeneous metapopulation networks," 
Phys. Rev. Lett., 99, 148701 (2007). 

[25] N M Ferguson, D A Cummings, C Fraser, et a/., "Strate- 
gies for containing an emerging influenza pandemic in 
Southeast Asia," Nature (London), 437, 209 (2005). 

[26] We assume here that individuals return to their base loca- 
tion k before travelling somewhere else. The initial model 
defined by Eqs. ^ allows to explore more general mo- 
bility patterns by relaxing this assumption. However, as 
long as the individual rates cj^^ are dominated by the 
dynamics that incorporate a base location, propagation 
speeds do not deviate significantly from the model we 
investigate here. 

[27] U Naether, E B Postnikov, and I M Sokolov, "Infection 
fronts in contact disease spread," Eur. Phys. J. B, 65, 
353 (2008). 

[28] F Ball, D Mollison, and G Scalia-Tomba, "Epidemics 
with two levels of mixing," Ann. Appl. Prob., 7, 46 
(1997). 



